program mkhis ! -- makes histgroms from list-mode data file.
  !
  ! DESCRIPTION:
  !
  ! There are two sisters of this: mkhisq and mkhisqp. They all
  ! produces the same results, i hope. You would prefer mkhisq if it's
  ! true. If it's not, you should believe the results of mkhis because
  ! it has the most redundant checking code for discarding invalid
  ! list-mode data.
  !
  !
  ! LICENSE AGREEMENT:
  !
  ! You agree that we don't agree to anything.
  !
  !
  ! DISCLAIMER OF WARRANTIES:
  !
  ! This software is provided to you "AS IS," and the authors and
  ! contributors disclaim any warranty or liability obligations to you
  ! of any kind, implied or express, including but not limited to the
  ! warranties of merchantability, fitness for a particular purpose,
  ! and non-infringement.
  ! 
  use suty
  use kuty
  implicit none
  interface
     subroutine decd(pe,pd,tad1,tad2)
       integer*8,intent(in),value::pe
       integer*4,intent(out)::pd
       integer*4,intent(out)::tad1(3)
       integer*4,intent(out)::tad2(3)
     end subroutine decd
  end interface
  !
  integer,parameter::nd=8         ! total number of blocks
  integer,parameter::ndh=nd/2     ! number of blocks in a plane
  integer,parameter::mgap=4       ! number of crystals in a gap
  integer,parameter::nt=64        ! number of crystals in transaxial
  integer,parameter::na=20        ! number of crystals in axial
  integer,parameter::nx=na*ndh+mgap*(ndh-1) ! number of crystals in x
  integer,parameter::ny=nt        ! number of crystals in y
  integer,parameter::MAXPATHLEN=1024
  !
  integer*8,parameter::fbs=1024*1024*1024 ! maxmum byte to read at once
  integer*4,allocatable::buf(:)
  integer*2,allocatable,target::histp(:,:,:,:)
  integer*2,allocatable,target::histd(:,:,:,:)
  integer*2,allocatable::bodyhist(:,:,:,:)
  integer*2,pointer::hist(:,:,:,:)
  character*(MAXPATHLEN) fin,str
  character*(MAXPATHLEN),target::fout,foutp,foutd
  character*(MAXPATHLEN),pointer::s
  integer istat,len
  integer*8 n,ie,b1,b2,fs
  integer*4,target::tad(3,2)
  integer*4,pointer::t1,a1,d1,t2,a2,d2
  integer*4 pd,tad_max(3,2),tad_min(3,2)
  integer if,nf,id,ia,narg  
  integer off(ndh)
  integer flg
  integer ix1,iy1,ix2,iy2
  logical lo_found
  integer*8 nb
  integer nok,npair,ndelay,nlost,ninv
  integer noverflow,nflip,nunderflow
  integer bphist(4,4)

  integer nee,neb,nbb

  narg=command_argument_count()
  if(narg<3) then
     call usage
     stop "*** too few arguments."
  end if
  
  fout=""
  foutp=""
  foutd=""
  nf=0
  nullify(s)

  do ia=narg-1,1,-2
     call get_command_argument(ia,str,len)
     if(len==2) then
        select case(str(2:2))
        case("o")
           s=>fout
        case("p")
           s=>foutp
        case("d")
           s=>foutd
        case default
           call usage
           stop "*** invalid option."
        end select
        if(associated(s)) then
           call get_command_argument(ia+1,s)
           nullify(s)
        end if
        nf=ia-1
     else
        if(fout/="".or.foutp/="".or.foutd/="") then
           exit
        else
           call usage
           stop "*** please specify outfile correctly."           
        end if
     end if
  end do

  tad=0

  do id=1,ndh
     off(id)=mgap*(id-1)
  end do

  allocate(histp(nx,ny,nx,ny))
  allocate(histd(nx,ny,nx,ny))
  allocate(bodyhist(nx,ny,nx,ny))
  histp=0
  histd=0
  
  ninv=0
  nok=0
  ndelay=0
  npair=0
  nlost=0
  nb=0
  noverflow=0
  nunderflow=0
  nflip=0
  t1=>tad(1,1); a1=>tad(2,1); d1=>tad(3,1)
  t2=>tad(1,2); a2=>tad(2,2); d2=>tad(3,2)              

  bphist=0

  nee=0
  neb=0
  nbb=0

  do if=1,nf

     call get_command_argument(if,fin,len)
     fin=trim(adjustl(fin))     
     write(*,"(a,x,$)") "Opening file: "//trim(fin)
     open(unit=1,file=fin,status="old",access="stream",iostat=istat)
     if(istat/=0) stop "*** cannot open file."
     fs=fsize(fin)
     write(*,*) trim(i2a(fs))//" byte"
     if(fs==0) then
        close(1)
        cycle
     end if

     n=fs/4
     if(allocated(buf)) deallocate(buf)
     allocate(buf(n),stat=istat)
     if(istat/=0) stop "*** cannot allocate buf"

     b1=1
     do ie=1,n/fbs+1
        b2=min(b1+fbs,n)
        read(1,iostat=istat) buf(b1:b2)    
        if(istat/=0) stop "*** error during read."
        nb=nb+(b2-b1+1)
        b1=b2+1
     end do
     close(1)

     lo_found=.false.
     do ie=1,n
        if(lo_found) then
           flg=get_flg(buf(ie))
           if(flg==0) then

              ! pair found
              lo_found=.false.
              npair=npair+1
              call decd(loc(buf(ie-1)),pd,tad(1,1),tad(1,2))

              tad_max=max(tad,tad_max)
              tad_min=min(tad,tad_min)

              !
              if(t1>=nt.or.t2>=nt.or.a1>=na.or.a2>=na) then
                 write(*,*) "*** WARNING: index overflow detected."
                 write(*,*) "t1,a1,t2,a2=",t1,a1,t2,a2
                 noverflow=noverflow+1
                 cycle
              end if
              if(t1<0.or.t2<0.or.a1<0.or.a2<0) then
                 write(*,*) "*** WARNING: index underflow detected."
                 write(*,*) "t1,a1,t2,a2=",t1,a1,t2,a2
                 nunderflow=nunderflow+1
                 cycle
              end if
              !
              if(pd==1) then
                 hist=>histp
              else
                 hist=>histd
              end if
              tad=tad+1 ! <<<

              if(d1<=ndh) then
                 if(d2<=ndh) then
                    ninv=ninv+1
                    cycle
                 end if
              else
                 if(d2>=ndh+1) then
                    ninv=ninv+1
                    cycle
                 end if
                 nflip=nflip+1
                 write(*,*) "*** WARNING: index flip detected."
                 write(*,*) "d1,d2=",d1,d2
                 CYCLE ! <<< dont know what to do but for safty
!!$                 t1=>tad(1,2); a1=>tad(2,2); d1=>tad(3,2)
!!$                 t2=>tad(1,1); a2=>tad(2,1); d2=>tad(3,1)
              end if

              if(d1>ndh.or.d2>nd) then
                 write(*,*) "*** WARNING: index overflow detected."
                 write(*,*) "d1,d2=",d1,d2
                 noverflow=noverflow+1
              end if
              if(d1<=0.or.d2<=ndh) then
                 write(*,*) "*** WARNING: index underflow detected."
                 write(*,*) "d1,d2=",d1,d2
                 nunderflow=nunderflow+1
              end if
              !
              ix1=(d1-1)*na+(na-a1+1)+off(d1)
              iy1=t1
              ix2=(d2-(ndh+1))*na+a2+off(d2-ndh)       
              iy2=t2
              if(ix1>nx.or.ix2>nx.or.iy1>ny.or.iy2>ny) then
                 write(*,*) "*** WARNING: index overflow detected."
                 write(*,*) "t1,a1,d1,t2,a2,d2=",t1,a1,d1,t2,a2,d2
                 write(*,*) "ix1,iy1,ix2,iy2=",ix1,iy1,ix2,iy2
                 noverflow=noverflow+1
                 cycle
              end if
              if(ix1<=0.or.ix2<=0.or.iy1<=0.or.iy2<=0) then
                 write(*,*) "*** WARNING: index underflow detected."
                 write(*,*) "ix1,iy1,ix2,iy2=",ix1,iy1,ix2,iy2
                 nunderflow=nunderflow+1
                 cycle
              end if
              hist(ix1,iy1,ix2,iy2)=hist(ix1,iy1,ix2,iy2)+1
              if(pd/=1) ndelay=ndelay+1
              nok=nok+1
              if(isEdge(a1,t1)) then
                 if(isEdge(a2,t2)) then
                    nee=nee+1
                 else
                    neb=neb+1
                 end if
              else
                 if(isEdge(a2,t2)) then
                    neb=neb+1
                 else
                    nbb=nbb+1
                    bodyhist(ix1,iy1,ix2,iy2)=bodyhist(ix1,iy1,ix2,iy2)+1
                 end if
              end if

              bphist(d1,d2-4)=bphist(d1,d2-4)+1

           else
              nlost=nlost+1
              if(flg==1) then
                 lo_found=.true.
              else
                 lo_found=.false.
              end if
           end if
        else        
           if(get_flg(buf(ie))==1) then
              lo_found=.true.
           else
              nlost=nlost+1
           end if
        end if
     end do
  end do

  write(*,*) "Number of files read:      ",nf
  write(*,*) "Number of bytes read:      ",nb*4
  write(*,*) "Number of total events:    ",npair
  write(*,*) "Number of overflow events: ",noverflow
  write(*,*) "Number of underflow events:",nunderflow
  write(*,*) "Number of fliped events:   ",nflip
  write(*,*) "Number of accept events:   ",nok
  write(*,*) "Number of random events:   ",ndelay
  write(*,*) "Number of invalid events:  ",ninv
  write(*,*) "Number of lost 4-bytes:    ",nlost

write(*,*) "Number of edge-edge events:",nee
write(*,*) "Number of body-body events:",nbb
write(*,*) "Number of edge-body events:",neb
 
  if(foutp/="") call write_hist(histp,foutp)
  if(foutd/="") call write_hist(histd,foutd)
  
  write(*,*) "Non-edge pair Max Prompt: ",maxval(bodyhist)
  write(*,*) "Max Prompt: ",maxval(histp)
  write(*,*) "Max Delayed:",maxval(histd)
  hist=>histp
  hist=histp-histd
  where(histp<0) histp=0
  write(*,*) "Max True:   ",maxval(hist)

  if(fout/="") call write_hist(hist,fout)
  
contains
  
  subroutine usage
    write(*,100)
100 format(x, "Usage: mkhis infiles",/,&
         8x,"{-o outfile_true | -p outfile_prompt | -d outfile_delayed}")
  end subroutine usage
  
  integer function get_flg(e)
    integer*4,intent(in)::e
    get_flg=ishft(e,-30)
  end function get_flg

  subroutine write_hist(h,outfile)
    integer*2,intent(in)::h(:,:,:,:)
    character*(*),intent(in)::outfile
    integer istat
    character*128 s
    write(*,"(a,$)") "Writing "//trim(outfile)
    open(unit=1,file=outfile,access="stream",status="replace",iostat=istat)
    if(istat/=0) then
       write(*,*) "*** error during open"
    else
       write(1,iostat=istat) h
       close(1)
       if(istat/=0) then
          write(*,*) "*** error during write"
       else
          s=trim(i2a(nx))//"x"//trim(i2a(ny))
          write(*,*) "size: "//trim(s)//"x"//trim(s)//"x"//trim(i2a(sizeof(h(1,1,1,1))))//" = "&
               //trim(i2a(sizeof(h)))//" bytes"
       end if
    end if
  end subroutine write_hist

  logical function isEdge(a,t)
    integer*4,intent(in)::a,t
    isEdge=((a==1.or.a==na).or.(t==1.or.t==nt))
  end function isEdge

end program mkhis
